function [model,errorMes] = solveATSM(params,setup)

%% Allocating memory 
maxMat     = max(setup.matSelect);            % max maturity
ny         = size(setup.matSelect,2);         % number of observables
nx         = setup.nx;                        % number of factors

% Unfold parameters in params to get the QTSM
if nx == 1
    alpha   = params.alpha;
	beta    = [params.beta1]';
    phi     = [params.phi11];
    mu      = [params.mu1     ]';
    sigma   = [params.sigma11 ];
elseif nx == 2
    alpha   = params.alpha;
	beta    = [params.beta1  params.beta2  ]';
    phi     = [params.phi11  params.phi12;
               params.phi21  params.phi22  ];
    mu      = [params.mu1    params.mu2    ]';
    sigma   = [params.sigma11  0                 ;
               params.sigma21  params.sigma22    ];
elseif nx == 3
    alpha   = params.alpha;
    beta    = [params.beta1  params.beta2   params.beta3]';
    phi     = [params.phi11   params.phi12   params.phi13;
               params.phi21   params.phi22   params.phi23;
               params.phi31   params.phi32   params.phi33];
    mu      = [params.mu1     params.mu2     params.mu3]';
    sigma   = [params.sigma11  0                 0              ;
               params.sigma21  params.sigma22    0              ;
               params.sigma31  params.sigma32    params.sigma33];
elseif nx == 4
    alpha   = params.alpha;
    beta    = [params.beta1  params.beta2   params.beta3   params.beta4]';
    phi     = [params.phi11  params.phi12   params.phi13   params.phi14;
               params.phi21  params.phi22   params.phi23   params.phi24;
               params.phi31  params.phi32   params.phi33   params.phi34;
               params.phi41  params.phi42   params.phi43   params.phi44];
    mu      = [params.mu1    params.mu2     params.mu3     params.mu4]';
    sigma   = [params.sigma11       0                  0                0;
               params.sigma21       params.sigma22     0                0;
               params.sigma31       params.sigma32     params.sigma33   0;
               params.sigma41       params.sigma42     params.sigma43   params.sigma44];
end

% Law of motion for factors under Q 
h0Q = phi*mu;
hxQ = eye(nx)-phi;

% Law of motion for factors under P
%h0  = h0Q + f0;     %note MPR is scaled by sigma
%hx  = hxQ + fx;

% Check if co-variance matrix for innovations is positive semi-definite
sigma2 = sigma*sigma';
if any(eig(sigma2)<=0)      
    errorMes = 1;
    model = [];
    return
end
%% Compute bond pricing coefficients
% Initialise memory for bond price coefficients in state 1
A = zeros(1,maxMat);
B = zeros(nx,maxMat);

% Compute recursive coefficients for bond prices
sigma2 = sigma*sigma';
for k = 1:maxMat   
    % One-period bond prices use the initial conditions
    if k == 1
        A(1,k) = -alpha;
        B(:,k) = -beta;

    % Longer bond prices depend on the previous maturity
    else
        A(1,k) = -alpha + A(1,k-1) + B(:,k-1)'*phi*mu + 0.5*B(:,k-1)'*sigma2*B(:,k-1);
        B(:,k) = -beta + (eye(nx)-phi)'*B(:,k-1); 
    end
end

%% Reporting output in the desired format: annualized yields in pct.
% Spot: y_{n,t} = -(1/n)*p_{n,t}
g0All= 12*(-A(1,:)'./(1:maxMat)');
gxAll= 12*(-B(:,:)'./repmat((1:maxMat)',1,nx));
matAll=1:maxMat;
g0   = 12*(-A(1,setup.matSelect)'./setup.matSelect');
gx   = 12*(-B(:,setup.matSelect)'./repmat(setup.matSelect',1,nx));
r0   = 12*alpha;
rx   = 12*beta';
% saving output
model = struct('g0',g0,'gx',gx,...
    'alpha',alpha,'beta',beta,'phi',phi,'h0Q',h0Q,'hxQ',hxQ,'sigma',sigma,'matSelect',setup.matSelect,...
    'r0',r0,'rx',rx,'g0All',g0All,'gxAll',gxAll,'matAll',matAll);

% Report no error if applicable
errorMes = 0;


end